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We develop a methodology for the calculation of surface free energies based on the probability 
distribution of a wandering interface. Using a simple extension of the NpT sampling, we allow 
the interface area to randomly probe the available space and evaluate the surface free energy from 
histogram analysis and the corresponding average. The method is suitable for studying systems with 
either continuous or discontinuous potentials, as it does not require explicit evaluation of the virial. 
The proposed algorithm is compared with known results for the surface tension of Lennard-Jones 
and Square Well fluid, as well as for the interface tension of a bead-spring polymer model and good 
agreement is found. We also calculate interfacial tensions of freely jointed tangent hard sphere chains 
on athermal walls for a wide range of chain lengths and densities. The results are compared with 
three different theoretical approaches, Scaled Particle Theory, the Yu and Wu density functional 
theory and an analytical approximation based on the latter approach. Whereas SPT only yields 
qualitative results, the last two approaches are found to yield very good agreement with simulations. 

PACS numbers: 68.03.Cd, 68.08.-p, 68.03.-g, 68.35.Md, 68.47.Mn 

INTRODUCTION 

Interfacial phenomena have been a matter of great research interest for centuries [l[ . With the growing importance 
of nanotechnology, however, this field of research is expected to become even more relevant. As the surface to volume 
ratios diminish in miniaturized devices, surface effects become increasingly important and strongly influence the state 
of the systems Q , including new and surprising behavior [3], 0] . 

Whereas the structure of the interface can be quite complex at the microscopic level 0, 0, 0, HI, @, EH II | > from 



a thermodynamic point of view only the interfacial tensions are really required to describe the system's behavior in 
most instances (e.g., Young's equation). Statistical mechanics provides a link for the calculation of interface tensions 
in terms of the average of a mechanical property, namely, the anisotropy of the pressure tensor fl2| . This route 
has been applied to estimate interface tensions, either theoretically [lj, or by means of computer simulations [l3j |. 
Results are known for some prototype systems, inclu ding the Lennard-Jones Square Well and Gay-Berne 
[l6| fluids, as well as bead-spring polymer chains [l7l Il8l Hil . Unfortunately, the explicit evaluation of the pressure 
tensor is already a very subtle matter for simple fluids [ll. Il7j. The presence of discontinuities in the potential, either 
inherent in the model or resulting from truncation in computer simulations can cause difficulties [lij ] , and at any rate 
need the approximate evaluation of Dirac delta functions. For molecular fluids the difficulties may become even more 
important, since the bonding is often rigid and the resulting contribution of the constraining forces becomes difficult 
to estimate [l9j . 

Several methods have been proposed in the literature to overcome the difficulties related to explicit evaluation of 
the virial. A useful methodology which is particularly interesting in the vicinity of the critical point was proposed long 
time ago by Binder [2(| , and has received renewed interest recently 0, [2l|, [22] • Although cast in terms of probability 
distributions, this method is equivalent to a full thermodynamic integration from the one phase state to a two phase 
state 23]. In that sense, it can be considered as an extension of thermodynamic integration employed to calculate 
interfacial tensions from adsorption isotherms (24|. The need for sampling all states between the homogeneous bulk 
phase and the phase coexisting states (slab geometry) makes it fairly time consuming away from the critical point. On 
the other hand, one gets a large amount of information on the intermediate metastable states which can be exploited 



for use in nucleation studies 25, 26]. Alternatively, one can explicitely simulate two phase states only and attempt 
to collect averages containing information on the interface. One such method is based on the capillary Hamiltonian 
approximation, whereby the spectra of interfacial capillary waves is analyzed, and the surface tension is obtained 
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from the long wavelength behavior @, [H, E3. This method may be applied only for the liquid-vapor interface and 
presents difficulties related to the explicit location of the interface boundary but does not seem to be affected by 
this arbitrariness at long wavelengths [28| . A third route that allows to calculate interfacial tensions without explicit 
evaluation of the virial is based on the perturbation of the two phase system, by increasing or decreasing the interface 
a finite amount. For curved interfaces, this method resembles the philosophy of Scaled Particle Theory 29|, 3(J 31 1, 



and was exploited by Bresme and Quirke in order to calculate the interfacial tension of colloidal nano particulates 
[32I l33l ]. Another perturbative approach, the Test Area Method of Gloor et al., has been recently applied to the 



study of flat interfaces, providing reliable estimates for spherical and anisotropic fluids |34|, |35|. These methods are 
special cases of a more general methodology which can be applied to measure first derivatives of the free energy [36| , 
with the Widom Test particle method as the most widely known application [37I ]. In the perturbative methods, the 
difficulty is the choice of surface increment, and the reliable average of the resulting small energy change, as is the 
case in the evaluation of numerical derivatives. In this work we propose a new method which has features of both 
the Binder method, in the sense that it relies on histogram analysis, and the perturbative approaches employed by 
Bresme and Quirke and Gloor et al., in the sense that the system is subject to small perturbations of the surface 
area. The difference is that the perturbations add up and the interface is allowed to wander freely. An analysis of the 
resulting probability distribution, or, alternatively, the corresponding average, yields the required interface tensions. 

Another route for calculating the interfacial tension is the density functional theory. In this approach the grand 
potential of the system is a functional of the local density, therefore the interfacial tension is readily accessible once the 
functional is minimized. Several density functionals and related self consistent field theories for polymeric fluids have 
been proposed in the literature [H, [H, [3!^, H(| ■ One class of functionals is based on the polymer reference interaction 
site model. 13, H, 41, 42 1 Another class is based on the Woodward formalism [1, |43 , [44 1 . In this work we use the 



functional belonging to the latter class and proposed by Yu and Wu [43} . This theory is based on Wertheim's first-order 
thermodynamic perturbations theory :45j (TPT1), and is one of the most accurate functionals for tangent hard sphere 
systems. However, besides the comparison with our new simulation methodology our choice of the functional serves 
another purpose. The Yu and Wu theory is also used as a starting point for derivation of an analytical approximation 
for the interface tension of tangent hard-sphere chains at hard walls. This simple approximation will turn out to be 
in a good agreement with the simulation data. 

The content of the paper is organized as follows. In next section we briefly consider the surface thermodynamics of 
systems with inhomogeneous density profiles along one direction perpendicular to the interface. This applies for two 
cases of interest, free liquid-vapor interfaces and fluid-substrate interfaces. Having considered this formal aspects, 
we present the new method proposed for the calculation of interface tensions in section III. Results for the different 
systems studied, from the Lennard-Jones fluid to athermal chains is presented in section IV, where comparison is 
made to DFT results and Scaled Particle Theory. Finally section V presents our conclusions. 

SURFACE THERMODYNAMICS OF SLIT PORE GEOMETRY 

In this section we will study the thermodynamics of systems having an inhomogeneous interface profile along a 
single direction, z. This includes two systems of interest, among others. One consists of (say) a liquid slab surrounded 
by vapor (or any other coexisting phase) in a system with periodic boundary conditions in all three directions. The 
other is a fluid in contact with a substrate, with periodic boundary conditions only in the directions perpendicular 
to z. For the sake of generality let the container have orthogonal shape, with sides of length L x , L y and L z (note 
that in a slit pore the choice of L z amounts to an arbitrary definition of the dividing surface). Work can be done 
on the system by changing the shape in either of two ways. Firstly, one can perform displacements of the xy plane 
along the perpendicular direction. Let p± be the external pressure exerted on this moving plane. The work done by 
an infinitesimal increment of L z is then given by dw = —L x L y p± dL z . Alternatively, one can perform perpendicular 
displacements of the yz and zx planes, thereby producing an infinitesimal deformation of the xy plane. The resulting 
amount of work is then dw = —L z p\\ d(L x L y ), where p|| is the average external pressure exerted on the yz or zx 
planes: 

1 f 1 

P " = L~ J 2 +Pyy( z )) dz (!) 

Accordingly we may write for the infinitesimal changes in energy: 

dU = T dS - A Pl _ dL z - L z p\\ dA + fx dN (2) 
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where A — L x L y . Introducing a grand potential by means of the Legendre transform ft = U — TS — fiN yields: 

dfl = -SdT - Nd^i- A P± dL z - L zP[l dA (3) 

In practice, we will be concerned with deformations which preserve the volume, hence AdL z + L z dA = 0. It then 
follows that 

In the limit of large L z , where the two interfaces do not interact any longer, the above derivative is identified with 
twice the interface tension, 700. For the more general case where interaction among the interfaces is possible, we may 
therefore write: 



dfi = —SdT — iVd;U + 2 r y(L z ) dA (5) 

Note that f2 is a first order homogeneous function of the interface area (cf. Eq. (|3J)). This allows us to write: 

fl = 2~/A-p x AL z (6) 

This equation will be exploited later on in order to calculate 7. 

In order to make contact with the standard thermodynamics of slit pores, two new definitions must be introduced. 
Firstly, a solvation force is defined as the excess of p± over the expected bulk pressure f s (L z ) = p±_ — p. Second, we 
introduce a surface free energy per unit area, 2u(L z ) — — (pii — p)L z . With these definitions, Eq. ([3]) then becomes: 

dfl = -SdT -pdV - Ndfi + 2ujdA- Af s dL z (7) 

The resemblance with the thermodynamics of slit pores is now apparent, although the volume here is just that of 
the actual system, V — L x L y L z with no mention of the reservoir. The equivalence of the two descriptions follows 
by considering the surface excess free energy. The free energy increment of a bulk system of equal shape and size is 
readily obtained either from Eq. ((3|) after replacement of N, S and pii and p± with the corresponding bulk values or 
from Eq. |(7J) by retaining the first three terms in the right hand side. In either case, we obtain the excess free energy 
differential that is standard in slit pore thermodynamics [46| . 

dAQ ex = -AS ex dT - AN ex d^i + 2uj dA - Af s dL z (8) 

where subscript ex stands for excess amounts over the corresponding bulk property. 

Whether one uses Eq. ([3]) or Eq. |(7J) as the fundamental relation is just a matter of convenience. In experimental 
work on slit pores, one has control over L z , A is difficult to modify and Eq. ([7]) is more convenient. In the case of the 
simulations performed in this work Eq. ([3]) and Eq. ([5]) will prove more useful since only the independent variables 
are made explicit and none of these refer to the reservoir, which cannot be actually simulated. Furthermore, we 
expect j(L z ) to converges fast. Indeed, to a good approximation the difference between j(L z ) and 700 arises from 
contributions well inside the infinitely large slit were pu — p± is zero anyway. 



NEW METHOD FOR THE CALCULATION OF INTERFACIAL TENSIONS 



Formal aspects 



The method proposed is conce ptu ally very simple, and resembles other known methods where the simulation's box 
shape is allowed to fluctuate 47, 4^, 49|, 5(J. In the absence of any other constraint, the shape of the system at fixed 



temperature, chemical potential and volume will result from a competition between p\\ and p±. The rules of this 
competition are dictated by a probability density, f^vT{A), and the corresponding free energy, as follows: 

^VT{A) = -k B T\nf^ VT {A) (9) 

This provides the link between the statistical mechanical property, / and the macroscopic observables. By use of 
Eq. ([6]), we find right away: 

Uvt(A) = e -P(^A- P± v) (10) 
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This is the key result of the method proposed: sampling the interface at constant volume from the above distribution 
will readily yield the interfacial tension. For systems where the interfaces do not interact and both 7 and p± are 
system size independent, it is clear that the slope of InfuVT yields the interfacial tension right away. For finite 
systems, f u VT the analysis is complicated by the non-trivial dependence of j{L z ) and also p±(L z ). If sufficiently fine 
data ara available, however, one can still recover j(L z ), since, according to Eq. ([5]) and Eq. ([5]), we have 

KhL-w*.) (") 

In practice, however, Eq. (fl"T)|) as stated above is hardly useful, because f^vri-A) is a monotonous function; i.e., the 
unconstrained system does not have an equilibrium state and during a simulation it will either stretch infinitely in 
order to eliminate the interface if 7 is positive, or increase its lateral size continuously if 7 is positive. 

In order to overcome this difficulty, two different solutions seem possible. One is to add an extra work contribution 
that counterbalances the interfacial tension (balance sampling), yielding a proper equilibrium state. The other one 
is to constraint the system's area within a bracketed interval chosen conveniently (bracketed sampling). In practice, 
both physical situations amount to the introduction of a modified or constrained probability density of the form: 

Uvt(A) = e-WWfuM*) (12) 

where W(A) is a suitably chosen function of A. 

This distribution may be obtained from knowledge of the molecular interactions by sampling a modified grand 
canonical probability density of the form: 

UVT(N, A; {V} N ) OC e ^N e -W{{r}N)+W { A)) (13) 

where U is the intermolecular potential energy, and {r}jy is the set of position vectors that define the microstate 
of a system with overall N molecules. An obvious choice for the evaluation of the above probability density is the 
standard Grand Canonical Monte Carlo Method, with the addition of an extra box shape sampling. The box shape 
sampling is performed in the same spirit as in the standard NpT simulations. Each position vector r i = (r x ,r y ,r z ) is 
transformed into a dimensionless position vector t i = (r x /L x ,r y /L y ,r z /L z ) and the relevant density distribution to 
sample the shape of the simulation box becomes: 

UVT{N, A; {t} N ) OC e ^N v N e -0(U({t} N -L x ,L y ,L^+W(A)) (14) 

where we have made explicit the parametric dependence of the intermolecular energy on the box shape that results 
from the transformation. 



Practical implementation 

In this section we discuss the implementation of the proposed method using Monte Carlo simulations for the 



sampling of Eq. (fl3|) . Other possible choices such as the Wang-Landau sampling 5l| , and Molecular Dynamics will 
not be discussed here but are also possible. 

A simulation box of suitable size is chosen with initial values of L z = D about twice those of L x = L y = L. 
The simulation is then carried out using a standard grand canonical method, organized in cycles of several regular 
canonical and grand canonical trials. After the end of each cycle, a box deformation is attempted as follows. First, an 
area increment A is chosen uniformly from the interval [—A max , A max ] and a new trial area A n = A + A is proposed. 

1 /2 

The new box parameters are set accordingly as L n — An and D n — V/A n - The attempted move is then accepted 
with probability: 

P accept = min ( 1, e -m,-u 0+ w n -w o) \ (15) 



Obviously, the choice of W is a crucial issue for the performance of the method. At first sight, a work contribution 
that competes against the ^A term so as to produce a well defined equilibrium state seems the most appropriate. 
Unfortunately, due to the small systems considered the equilibrium value of A must fall within a fairly limited interval. 
Whereas this is can be achieved in principle, in practice it amounts to a good a priori knowledge of 7. Therefore, after 
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several attempts with possible W functions, it was found that the most convenient procedure is simply to bracket the 
interface area within an interval chosen beforehand, by means of the following W function: 

A <c A m i n 

E -PW{A) = I 1 <A< (lg) 

A ^> A max 

Within the chosen interval, the method produces the probability distribution of Eq. (|10p . For non-interacting inter- 
faces, a logarithmic plot of J^vt{A) yields the interfacial tension right away, while for interacting interfaces, Eq. (fTTJ) 
would have to be employed in order to determine j(L z ). In the former case, routine calculation of interfacial tensions 
may become somewhat cumbersome, and a simpler data analysis might be desirable. One such possibility is to exploit 
the expectation value of the surface area, which reads: 

(A) A 1/2 1 1 1 + e- 2 ^ 

AA 2pjAA 21-e- 2 ^ AA [ ' 

where A A = A max — A m i n and Ai = 1 /2{A max + A m i n ). The root of the above equation is the maximum likelihood 
estimate of 7 given (A) and may be solved by means of a Newton-Raphson method. The first order solution, 
2(3jAA — \2{Ai/2 — (A))/AA, may be used as an excellent first guess, producing convergence in less than four 
iterations. In fact, this approximation is good in most practical situations, as it produces deviations that are less than 
3% from the exact solution for values of (A1/2 — (A))/AA in the range [—0.1,0.1]. 

So far, the methodology described applies for the calculation of interfacial tensions in an open system. The choice of 
variables, /i, V and T is appropriate for calculating tensions of a fluid against a wall. Interfacial tensions between two 
coexisting phases are also of great importance, but the thermodynamic description employed here is not applicable 
right away. An initial state with two such phases will spontaneously transform into one or the other phase and 
completely eliminate the interface. This problem is solved immediately simply by replacing /1 with N as the natural 
variable, with the number of particles fixed within the immiscibility gap. All the framework that was employed in 
this section is then immediately applicable by simply replacing the grand free energy by the Helmholtz free energy, 
and ignoring the sampling over the particle number. 

MODEL AND RESULTS 
Molecular models studied 

The bracketed interfacial area sampling described in the previous section will be applied to the study of four different 
systems. Firstly, we will consider an argon like model, with pair interactions that only depend on the mutual distance, 
r between the atoms and obey the usual Lennard-Jones potential: 

Secondly, a central force potential of the square well type will also be considered: 

00 r < a 

Vsw(r) = { -e a<r<Xa (19) 
r > Act 

where A defines the width of the square well and a and e are the usual range and energy parameters of the models. 

Thirdly, we will consider a well known bead-spring polymer model which has been extensively studied. In this 
model, all beads interact with each other, whether in the same or different polymer by means of a truncated and 
shifted Lennard-Jones potential with R c — 2 ■ 2 1 / 6 a\ Additionally, adjacent beads within the same polymer are held 
together by means of a Finite Extensible non Elastic (FENE) potential of the form: 

2 

$(r) = -fc.J&ln(l-4-) (2°) 

where k s = 15e/a 2 plays the role of spring constant and Roo — 1.5a is the maximum allowed displacement between 
bonded beads. Whereas for the simple atomic fluids only the vapor-liquid interfacial tensions will be calculated, for 
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the polymer we will also report wall-fluid interface tensions. The polymers will thus be confined between two parallel 
plates in the usual slit pore geometry, with fluid-substrate interactions of the form: 




where A plays the role of a Hamaker constant. 

The last model considered in this work will be a chain of freely jointed tangent hard sphere chains. In this case, 
beads interact with each other via a hard sphere potential of the form 

w»(o (22) 

Monomer within the same chain have a fixed bond length of a, and are thus bonded tangentially. The only other 
constraint is the non-local hard sphere interaction between beads more than one segment apart, so that the chain is 
otherwise fully flexible. We will consider systems with slit pore geometry and a purely repulsive and athermal wall of 
the form: 

t r / \ / oo z < a/2 

W*) = | 2 > CT/2 (23) 



Liquid— Vapor surface tensions 

Surface tension of the Lennard-J ones Fluid 

The methodology proposed was first explored and tested using the Lennard-Jones fluid, since several other methods 
have already been tested for this system. The interactions were truncated at R c = 2.5a, with no further shift of the 
potential (this corresponds to the spherically truncated ST Lennard-Jones model of Trokhymchuk and Alejandre 
[LJ|). The simulated boxes had an initial size of 13x13x45, containing a total of 2137 particles. The production stage 
was organized in cycles. Each such cycle consisted of 2000 attempted canonical attempts to displace a particle. Half 
of these attempts were standard metropolis displacements set to yield about 50% acceptance. The remaining half 
consist of biased deletion-insertion movements [521 ] . First, a particle chosen at random was deleted. The particle 
was then inserted in one among several randomly chosen positions with a probability given by the corresponding 
Rosenbluth weights. The attempted displacement was accepted or rejected following the configurational bias rules 
[53l ]. It is expected that this trial move will accelerate the proper equilibration between the vapor and liquid phases 
and provide a better sampling of the interface than the metropolis movements alone. After the end of each cycle, a 
box deformation attempt was performed, with A att set during the equilibration period to yield about 50% acceptance. 
Averages were collected over about half million cycles to one million cycles. 

Before considering the performance of the method proposed in the previous section, with the specific bracketing 
choice for W (cf. Eq. fl!6|) ). it is interesting to consider the performance of an alternative counter tension function 
selected so as to produce a proper equilibrium state. The surface tension will favor states such that jA is minimum. 
In order to equilibrate this contribution, an extra term acting in the opposite direction is needed. A possible choice 
that was considered is: 

W(A) = 2fA 2 /A (24) 

where Aq — V 2 I Z is introduced in the equation so that the parameter £ has the same dimensions as 7. At constant 
N, V and T, the relevant Helmholtz free energy for the system in this case is simply: 

F(A) = 2fA 2 /A + 2 1 A (25) 

Minimizing the free energy it is found that the equilibrium state of the surface area is such that: 



A_ _/^ 1/a 
A Q ' 



(26) 
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Therefore, the average interfacial area obtained from a simulation allows to calculate the surface tension. An obvious 
way to test the above equation is by performing several simulations with different £ values. A plot of (A/Aq) 2 
against £ will produce a straight line with slope 7. This test is shown in the inset of FigQ] for the ST-L J model at 
T = 0.85e/kB- The plot clearly shows a straight line with slope f3j — (0.779 ± 0.0I0)cr 2 , in good agreement with 
the /37 = (0.772 ± 0.010)er 2 result suggested by Trokhymchuk and Alejandre from explicit evaluation of the virial 
In that figure, results obtained by the balance method for several other temperatures is also shown. Good 
agreement with Refflil is found in most cases. A more detailed analysis of the results, together with the average 
surface area and £ values employed is reported in Table [B This is an important parameter, because 7 is known to 
show a significant dependence on the lateral system size [HH HE]. Our results were performed for similar system 
sizes as those reported in Refll4l and the comparison shows good agreement. The only large discrepancy is found 
for the relatively high temperature T = 1.127e/fes, close to the critical point. A reliable estimate of 7 at such high 
temperature is difficult for several reasons. The interface width diverges close to the critical point, the bulk densities 
are attained only asymptotically slow and the results will depend also on the perpendicular length L z and the length 
between the two interfaces as well. One advantage of the method proposed here is that the interface is allowed to 
fluctuate, so that the value of 7 that is obtained is averaged out over a finite interval of lateral lengths. For high 
temperatures these fluctuations can be large. Indeed, the interfacial areas sampled during the simulation of the system 
at T — f.l27e/fcs for £ = 0.0075e/<7 2 fall broadly in an interval of lateral lengths between 10 and 14o\ Unfortunately, 
these large fluctuations could actually be more a matter of concern than an advantage. For a very large system, the 
fluctuations about the equilibrium interface area will be Gaussian, and the average interfacial area will be equal to the 
extremum of the distribution. For finite systems as in our simulations, the large fluctuations imply that the interface 
area is actually sampling a non Gaussian skewed distribution of the form f(x) ~ exp(— (ax + b/x)), with x the random 
variable and a and b some coefficients. In such a case, the average and the extremum may deviate significantly, and 
Eq. ([25]) might not be a reliable estimate of the surface tension. 

However, the most important problem with the counter tension methodology is in fact the choice of £. Note that 
the average surface area that resulted from our simulations is such that the lateral length remains close to the initial 
value L = 13(7 for most temperatures (cf. Table HJ. Obviously, a look at Eq. (|25|) shows that this cannot be achieved 
accidentally. For an arbitrary choice of £, the resulting equilibrium geometry may produce problems. For any £ larger 
than 7, it is clear that the average equilibrium surface is larger than Aq. As a result, the lateral length becomes larger 
than the perpendicular length, and the system becomes metastable: i.e., the liquid slab will have tendency to rotate 
by 7r/2 so as to achieve a minimal surface area. For values of £ that are too large, on the contrary, the simulation box 
will stretch very much and produce a long slab with very small surface area. The surface tension that is obtained will 
then strongly depend on the equilibrium state, while use of Eq. (|25p which was derived by ignoring the system size 
dependence of 7 could be again not adequate. Obviously, one can always choose £ such that (A) falls within some 
desired range as it was done in this work, but this implies a reasonable first guess of 7, which is not always at hand. 

The bracketed interfacial area sampling proposed in the previous section overcomes all of these problems. Firstly, 
the sampling interval may be set at will right away, avoiding visits to undesired regions. Secondly, no prior knowledge 
of 7 is required. Table U shows the results obtained using the bracketing sampling strategy. Overall good agreement 
with both the virial and the balance sampling method is found, though the results for the near to critical temperature 
T = 1.127e/fcs are again in conflict with those obtained from the virial method. On the other hand, both methods 
employed in this work provide results in fair agreement for that temperature, given that the lateral system sizes 
studied differ considerably and the counter tension method could suffer from the problems discussed above. 

One important issue with the bracket sampling is the choice of bracketing interval [A m i n , A max \. On the one hand, 
a large interval is difficult to explore, because the states with large surface area are exponentially suppressed and so 
will suffer from poor sampling. On the other hand, large Ay and Ace values will lead to a small relative error of the 
slope, Ay/ Ax. Clearly, choosing the optimal value is not a trivial matter, although we expect the optimal choice will 
depend on the dimensionless parameter £ = 2(3jAA. For T = 0.70e/fcs and T = LOOe/Ais two set of simulations 
with different values of £ were performed. In set S, a small interval was chosen, corresponding to a small expected £ 
of about I. On the contrary, in set L we chose £ ~ 9 such that Jnvt(A) varied about 3 orders of magnitude. The 
statistical analysis based on partial results of 2/3j obtained from Eq. (jTTJ) after averaging the surface over 500,000 
cycles remain inconclusive. The error as measured from the standard deviation of the sample is very similar in 
both sets, and no significant differences in the mean are observed either. For the state point at T = 0.70e/fcs set 
S yields 2/57 = (2.245 ± 0.050)cr~ 2 while two independent runs of set L produced 2/3j = (2.282 ± 0.040)cr~ 2 and 
2/3 7 = (2.277 ± 0.027)cr~ 2 . Similarly, for T = 1 .00e/k B set S produced 2/3 7 = (0.490 ± 0.015)cr~ 2 while set L yield 
2/3 7 = (0.498 ± 0.017)cr- 2 and 2/3 7 = (0.490 ± O.Olljo-- 2 . 
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Surface tension of the Square— Well fluid 

Calculation of the surface tension from the anisotropy of the pressure tensor is particularly difficult for fluids with 
discontinuous potentials. In such cases, there appears a delta Dirac in the virial expression for every discontinuity 
in the potential. In practice, this means that the instantaneous virial can only be measured approximately and that 
the computer program must be tailored for each model. The method of interface sampling has the advantage of not 
requiring explicit evaluation of the virial, and this may prove particularly advantageous in this case. 

In this section the simple square-well fluid with A = 1.5a is chosen as a prototypical discontinuous model potential. 
Simulations were performed on systems containing overall 2222 atoms and initial system size of 13x13x45. The 
sampling strategy is exactly the same as with the Lennard-Jones fluid, with cycles of 2222 center of mass and 
configurational bias moves in the ratio 50:50. 

The difficulties in sampling the square well fluid interface are illustrated by considering the counter tension method 
explained in the previous section. In order to achieve a reasonable acceptance rate for the box deformations, the 
maximal random displacement has to be set to very small values, and it becomes difficult for the system to reach 
a meaningful equilibrium value of the interface area. For a temperature of T = LOOe/Zcs a value of A = 0.035cr 2 
was required to achieve the prescribed 50% acceptance ratio in box deformations. On the other hand, a value of 
A = 0.65fi 2 was obtained for the Lennard-Jones fluid at the lowest temperature studied. Assuming a simple random 
walk for the interface displacements, this would imply that the square well fluid needs (0.65/0.035) 2 more cycles than 
the Lennard-Jones in order to sample the same surface area. In this conditions, the counter tension method becomes 
inefficient and the bracket sampling strategy with a small interval of about 0.5<r 2 is more convenient. Even so, the 
simulations were performed over ten million cycles in order to achieve error bars comparable to those obtained for 
the Lennard-Jones. 

The results obtained using the bracket sampling arc shown in table [TTJ Good agreement with previous results using 



the virial method, the Test Area Method and Binder's method is found [15l. |34|. l56j. 

In the Test Area method, a simulation is carried out in an NVT ensemble of fixed shape. Every cycle, a 'test' 
deformation at constant volume is performed, and the corresponding change in configurational energy, AU is measured. 
Following Gloor et a!., it can be shown that the surface tension is given as: 

ln(e^ AC/ ) n 

^ = &o- kBT AA (27) 

where the subscript "0" denotes that the interface is sampled over the unperturbed system. In the limit where 
the method becomes exact, it becomes equivalent to measuring the anisotropy of the pressure tensor, with the virial 
evaluated numerically instead of analytically. Indeed, for small enough deformations, the instantaneous energy change 
may be written as: 

»im=f«) A, (28) 



V dA 



v 



Accordingly, an asymptotically exact Taylor expansion of Eq. (|27p yields right away 

f dU({r N }) 



7fa = (( q A " ) ) (29) 



vi o 



From the canonical partition function one can show that the analytical derivative in this equation is just the anisotropy 
of the pressure tensor. Both the TAM method and the interface sampling proposed in this work have the advantage 
of avoiding the explicit calculation of the virial. However, the interface sampling method could avoid difficulties 
encountered with the TAM method due to the nonequivalence between compressing and expanding the interface in 
discontinuous systems [34j |. For very small deformations the three methods should become equivalent, although the 
virial method would seem to yield results with smaller error bars [35| . 



Fluid— substrate interfacial tension of Lennard-Jones chains 



In this section the interfacial tensions of the bead-spring polymer model described at the begining of this section are 
considered. The interfacial properties of chains with 10 beads at a subcritical temperature of T = 1.68e/ks have been 
studied previously and are well known (il. Il7j|. The liquid-vapor surface tension has been calculated using Binder's 
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method, while the spreading coefficients ~/ wv — j w i were also evaluated with a related method all the way from the 
wetting to the drying transition . The interfacial tensions of the glassy polymer in narrow slit pores have also been 



measured using the virial method [17 1 



Simulations for the polymer model were carried out in the grand canonical ensemble. The chemical potential 
was fixed to its coexistence value and the initial box size was set to 13.8x13.8x27.6 (i.e., for the liquid state, this 
amounts to about 315 polymers in the simulation box). The configurational space was sampled by means of center of 
mass displacements, standard canonical configurational bias displacements and grand canonical configurational bias 
insertion/deletion attempts in the ratio 25:25:50. For the vapor phase only the last two kinds of trial moves were 
attempted. The maximum allowed box deformation attempts were adjusted during the equilibration period so as 
to provide roughly a 50% acceptance ratio. Averages were collected over roughly 1.5 million cycles for wall-liquid 
tensions and about 5 million cycles for the wall-vapor tensions. 

The results obtained for both the wall-liquid and wall-vapor interface tensions may be found in Table IIIII Also 
included are the differences j wv — j w i previously calculated using two different methods. The comparison shows excel- 
lent agreement for the spreading coefficient and thus provides strong support for the methodology. The liquid-vapor 
surface tension was also calculated using the interface sampling in the canonical ensemble with a 13.5x13.5x50 box 
geometry and 185 molecules. The result, 7;,, = (0.156±0.008)e/cr 2 is in good agreement with previous measurements 
using Binder's method, 7;„ = 0.160e/a 2 [9j. 

For many practical applications, the relevant property is not really the surface tensions alone, but the correspond- 
ing spreading coefficient, which governs both the wetting behavior of a single interface and the capillary condensa- 
tion/evaporation phenomenology. Figure [3] shows the results of the spreading coefficient as a function of e w . The 
points were the spreading coefficient cross the liquid-vapor surface tensions indicate the location of wetting/drying 
transitions. Independent measurement of the tensions remains interesting, because it allows to compare the free energy 
of the wall-fluid interface relative to that of a bulk system of equal volume. Such data is useful for comparison with 
self-consistent and density functional theories, but can be also relevant in process were the interface is free to grow, 
i.e., nucleation or self assembly of biological systems. Table Hill shows that both ^ wv and j w i are negative for strongly 
attractive walls but change sign as the effective Hamaker constant decreases. For positive but small values of e w the 
weak attractive interactions are not able to compensate for the free energy cost of a inhomogencous density profile. 
The fact that the liquid phase requires a wall strength larger than e w — 2e before it becomes negative is consistent 
with previously observed density profiles, which remain considerably depleted for e w as large as 3e and only start to 
develop monomer like oscillatory behavior at about that value. Actually, j w i is about one order of magnitude larger 
than j wv all the way from e w = to e w = 3.30, and dominates completely the behavior of ^ wv — j w i. At e w — both 
the wall-vapor and wall-liquid interfaces are unfavorable relative to the corresponding bulk phases, and the drying 
transition is actually driven by the very unfavorable wall-liquid interface. 



Interfacial tension of chains of tangent hard spheres at an athermal wall 

Having studied the performance of the method for fluids with discontinuous potentials (square well) and for polymer 
chains (FENE-LJ), we now consider a fluid made of tangent hard sphere chains adsorbed on purely repulsive athermal 
walls. This system has been subject of great interest in the last decade, as a p rot otyp ical model for the study of 
polymer-wall interfaces, as well as a test bead for density functional theories [24 1571. l58|. 

We have considered chains of 2, 4, 8, 12 and 16 monomers adsorbed on a hard wall. For each chain, we have 
calculated the interface tensions corresponding to six different bulk densities of approximately 0.1, 0.2, 0.3, 0.4, 0.5 
and 0.6 monomers per a 3 . We first performed bulk NVT simulations for those densities, and evaluated the chemical 
potential by means of Widom's test particle method. For the largest chains at heigh density, Widom's test failed 
and the grand potential could not be properly evaluated. An arbitrary value was imposed in those cases, based 
on the approximation that the excess chemical potential is linear in the chain length. The values of the chemical 
potentials employed in the simulations are collected in Tablc lTVl Once the desired chemical potential was known, grand 
canonical simulations were performed on a system with initial size of 15x15x35. The grand canonical distribution 
function was sampled by using center of mass, rotation, configurational bias and grand canonical configurational bias 
insertion/deletion moves in the ratio 15:15:35:35, except for dimers, were the ratio was set to 15:15:0:70. After suitable 
equilibration, averages were collected over 3 to 7 million cycles, with longer runs required for the highest densities. 

The choice of bracketing interval in this case is also a matter of concern, since the conditions change over a wide 
interval, and the trial displacements required to achieve 50% acceptance in the box shape moves varies from about 
3er 2 to about 3 x 10~ 2 <r 2 as the monomer density increases from 0.1er~ 3 to 0.6a~ 3 (this parameter being much less 
sensitive to chain length) . At this stage we developed a systematic method to decide the size of the lateral size sampling 



10 



interval which works well. A short simulation of about 100000 cycles is performed, and the lateral size mean squared 
displacement is calculated. A plot as a function of Monte Carlo cycles shows an initial linear and steady growth 
corresponding to the unconstrained diffusion process, followed by a plateau occasioned by the reflecting boundary 
conditions on the lateral size (Fig. QJ. Measuring an effective diffusion coefficient, D e ff, from the slope of the mean 
squared displacement within the first regime allows to estimate the expected lateral size displacement within a given 
number of cycles. Obviously, the lateral size needs to have reflected several many times if meaningful averages are to 
be collected. Therefore, the size of the sampling interval required for the condition to be obeyed within a prescribed 
amount, rtbiock of MC cycles may be estimated from D e ff. In practice, we chose the sampling interval from the 
distance traveled during 10000 cycles. 

The results of our calculations are presented in Table [V] As seen from the table, the surface tensions could be 
evaluated with relative errors lying between 1 and 10%, with accuracy deteriorating at the highest densities. Inspection 
of the results as plotted in Figf5] shows the interface tensions are always positive, increase with monomer density and 
decrease with chain length. Simulation results suggest an asymptotic behavior for long chains at low densities. For 
higher densities, however, the error bars are too large and this trend cannot be confirmed. 

Hooper et al. have calculated interface tensions of athermal chains with no intramolecular excluded volume inter- 
actions 24j . These authors find that the interface tensions of a 20-mer first increases with density, but then gradually 
decreases and actually becomes negative at heigh densities. The different behavior is related to the choice of dividing 
surface (cf. Eq. ([23|) . In this work the dividing surface coincides with the actuall wall surface (i.e., z = 0), while 
Hooper ct al. assume that the dividing surface is located at the distance of the closest approach of a segment to 
the surface (i.e., z — a/2). Our results can be easily recalculated to another choice of the dividing surface by using 
the simple relation 7 z=(T /2 = 7 2 =o — ^pa^wheie the subscripts stand for the location of the dividing surface and p is 



the bulk pressure (cf. Eq. (4.3) in Ref. |24T ). This equation shows that our choice of dividing surface yields what is 
sometimes known as the intermolecular contribution to the interfacial tension. For the proposed method, the choice 
of dividing surface is implied in the volume definition. [3^] Here we have employed V = A x L z , while restricting 
the system's volume to that space available to the molecule's center of mass we would require V = A x (L z — a) 
(a slit pore is used in practice, so \a must be substracted twice). These definitions dictate the length to breadth 
ratio of attempted deformations in the Monte Carlo simulation, and thus produce interfacial tensions according to the 
choice of dividing surface. In order to illustrate this point, we have performed some simulations for chains of length 
to = 12 and a dividing surface located at z = a/2. Fig{6] presents a plot of the simulation results, clearly showing 
the non-monotomic trend observed by Hooper et al • Included as squared symbols are results for j z=a /2 obtained 
from our previously calculated 7 z= o and shifted by — \pa. The results agree within the error interval. 

Along with the simulation results, we also include predictions from three different theoretical treatments. The first 
approach is based on Scaled Particle Theory [2§| , which has been shown to be accurate for fluids of spherical symmetry, 
whether hard spheres [35l . |59| , or Lennard-Jones j60| . In this approach one considers the free energy change that 
results from inserting a large par ticle, and identifies the chemical potential with macroscopic contributions arising 
from pV and surface work 



d/i ff = 2nja da + -^P& 2 da (30) 

where p a denotes the chemical potential of the large particle of diameter a inserted in the fluid. The macroscopic 
interface tension for a flat interface is obtained from the above equation in the limit of infinite dilution and particle 
diameter as a term proportional to a. Here we consider a fluid mixture of tangent hard sphere chains and large hard 
sphere particles as described by Wertheim's perturbation theorv[45j|. together with the Boublik equation of state for 
the hard sphere reference mixture [62| . After some tedious algebra, we obtain: 

2 1 f (3 - 6t7 2 + 37/ 4 ) ln(l - 7?) + 37/ + 14t? 2 - 5t7 4 - if + t? 6 ) 

nP<T 7 ~ (l + 77) 2 (l -r?) 2 1 +(37/- 57 7 2 + 27 ? 4 + 7; 5 -7 ? 6 )/m J 



The second approach is based on microscopic density functional theory for polymeric fluids [g, |43j, |44| . Within this 
framework the grand potential of the system is a functional of the local density of polymer p(R), 

fi[p(R)] = J dRp(R)(V ext (R) - M ) + /T 1 J dRp(R)(ln(p(R)) - 1) 

dRV b (R)p(R) + F ex [ PPS (r)} . (32) 



In the above R = (ri,r2, . . . , r m ) denotes the set of segment positions, V ex t is the external potential, while V& is 
the bonding potential that satisfies exp[— /3VJ,(R)] = I1™=l (^(l r i+i — r i| — a ))l '(47r<7 2 ). F ex [pps(r)] is the excess free 
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energy approximated as a functional of the average segment density denned as /9ps(r) = YliLi JdRS(r — rj)p(R). In 
this work we apply the Yu and Wu functional for F ex , which is based on fundamental measures' theory of Rosenfeld 
"43I . 6^ . According to this approach the excess free energy is a volume integral (3F ex = J dr&. The excess free energy 



density <f> is a simple function of a set of weighted densities {n a }, a = 3,2, 1, 0, VI, V2 defined as 

n a (r) = [dr'p PS (r')w a (r - r') . (33) 



The weight functions w a depend on geometrical properties of the segments and are given explicitely in Refs. 43, 6^. 

In the Yu and Wu theory the excess free energy density is represented as a sum of two contributions, $ = <&hs + ®p- 
$hs describes the reference mixture of hard spheres and is based on the Boublik equation of state 62. 64, 65| 



, , nxn 2 - n vl ■ n V2 
®hs = -n ln(l - n 3 ) H 

1 - m 

,/3 o n 3 + (1 - n 3 ) 2 ln(l - n 3 ) 

+ {nt, - 3n 2 n V 2 ■ n V2 ) — -. - 2 , , (34) 

367r(n 3 )^(l - n 3 ) 2 

while the excess free energy density due to the chain connectivity <I>p is an "inhomogencous counterpart" of the 
perturbation term in TPT1 

1 — 777 

$p = n (ln[y H s] , (35) 

m 

where £ = 1 — n\/ 2 ■ iiy 2 / (n 2 ) 2 . yns is the expression for the contact value of the hard-sphere radial distribution 
function 

1 n 2 aC, {n 2 <j) 2 C, 

l-n 3 4(1 -n 3 y 72(1 - n 3 ) 6 

The functional outlined above can be minimized numerically using the variational principle S£l/Sp(R) = 0. The 
resulting equilibrium density profile can be inserted into Eq. (|32[) yielding the grand potential which can be used to 
calculate the interfacial tension. 

The third theoretical approach utilizes the bulk limit of the microscopic density functional theory outlined above. 
Starting point is the same homogenous one-component polymeric fluid as in the SPT route, and again we consider 
the change of the grand potential, /3Af2, due to the insertion of a single big hard sphere into the system. This change 
can be expressed in terms of the derivatives of the excess free energy density 

9$ d$ 9$ 9$ 
(3AQ = P- Ca + ^- (2 + ^ Ci + jr- Co , (37) 
on 3 on 2 oni ouq 



where d i — 3,2,1,0 are characteristic functions of the shape of the inserted big particle 31, H, in]]- In particular £ 



,3 



and (2 are respectively the volume and the area of the big inserted hard sphere. If the functional was fully consistent 
with the SPT theory, the volume-dependent term d<fr/dn 3 should give the pressure, while the area-dependent term 
d&/dn 2 should give the planar wall- fluid interfacial tension (dQ/dni and dQ/dno should give the curvature-dependent 
corrections to the surface tension[3l|, 6(| 6?]])- The functional given by Eqs. (|34H36j) is not fully consistent with SPT 



but accepting this deficiency we can still identify the area-dependent term with the planar wall-fluid surface tension. 
Noting that in the bulk limit the vector weighted densities riv 2 and nyi vanish and n 3 — > 77, n 2 — > 6r]/a, n\ — > 3r]/a 2 ir, 
77 — * 6r]/a 3 ir the interfacial tension for tangent hard-sphere-polymer fluid at a hard wall is given as 

2 2 a$ 377(2 - 77) l-m2r7(3-77) 

7 = ™ = WW + 3 ( " " } + 777 (4-277) ' (38) 

The first two terms in the r.h.s of Eg 1381 represent the interface tension for the hard-sphere fluid at a hard wall and 
is the same as in the SPT approach presented above. The last contribution arises due to the chain connectivity but 
it differs from that of the SPT approach. Figure 5 presents the results for the interfacial tension of tangent hard 
sphere chains at a hard wall as a function of the chain length 777. The sets of curves were calculated for the bulk 
segment densities p\ = 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 (from bottom to top, respectively). The circles denote the 
interface sampling simulation data. The solid lines are the results obtained by numerical minimization of the DFT, 
while the dashed and dotted lines denote, respectively, the SPT approach (Eq. (|3T|) ) and the DFT-based analytical 
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approximation (Eq. (|38|) ). We note that the DFT results are in overall good agreement with simulations. The 
SPT approach always underestimates the interface tension and the agreement is only qualitative. The DFT-based 
analytical approximation gives poor results at low density but it significantly improves at higher density. Note that, 
as mentioned earlier, the interface tension is always positive because the dividing surface is set at the actual wall 
surface. This choice is particularly useful for the theoretical treatment, since adding the — \v<J extra contribution 
consistent with the choice of dividing surface at z = ier would obviously yield a less tractable equation. Figure [5] 
shows results for the interface tensions that result from this definition for chains of length m — 12. Consistent with 
our previous observations, all three theoretical treatments provide qualitative agreement with simulations, though the 
Yu and Wu DFT theory and the corresponding DFT-based analytical approximation are in better agreement than 
the SPT approach. 



CONCLUSIONS 



We have developed a method for the efficient calculation of interface tensions by means of computer simulations. 
The method allows the interfacial area to wander randomly and extracts the interface tension from the resulting 
probability distribution. This strategy shares common features with Binder's method and perturbative methods such 



as those proposed by Bresme and Quirke and Gloor et al [20|, |33j, [34[ . On the one hand it relies on histogram analysis 
or the related average of a sampled probability distribution; on the other hand the sampling is the result of small 
finite deformations of the simulation box. The method avoids the explicit calculation of the pressure tensor and is 
therefore of great generality. It can be applied likewise for simple or complex fluids with continuous or discontinuous 
potentials, and may be employed for interface tension calculations of free or bound interfaces. Any sampling method 
for the gand canonical or canonical distributions is appropriate, and improved sampling over wide ranges of box shape 
could be implemented using known techniques such as multicanonical (68j , successive sampling [6^ | , transition matrix 
[7(ij |. or Wang-Landau sampling (5l| . 

The proposed method was tested for the Lennard- Jones and Square Well fluids and good agreement was obtained 
with results found in the literature for similar system sizes. For systems with continuous potential, our results show 
that the method yields similar error bars at the same computational cost. The method runs into problems when 
the surface area cannot be sampled efficiently over a large interval. This is the case for the Square Well fluid, were 
somewhat greater error bars were found compared to the virial method at similar computational cost. On the other 
hand, the method is advantageous for the calculation of wall-substrate interfacial tensions at low density. In such 
cases, a large surface area can be sampled easily, and low error bars result. In such cases, calculation of the interface 
tension via the pressure tensor anisotropy is a problem, because the small number of collision events prevents accurate 
estimation of the virial. 

The technique was applied to the calculation of wall-fluid interfacial tensions for a bead-spring Lennard Jones 
chain, providing results for both the wall-vapor and wall-liquid tensions. The related spreading coefficient was found 
to be in good agreement with previous results Q . We also studied a system of tangent hard-sphere chains adsorbed 
on athermal walls for a wide range of chain lengths and densities. At constant monomer density the interface tension 
approaches an asymptotic constant value from above. For constant chain length, the interface tension increases 
regularly with density. The results where used to test three different theories. Density Functional Theory as proposed 
by Yu and Wu [43] , an analytical approximation based on that theory and another analytic result inspired on Scaled 
Particle Theory [291 ]. The latter approach reproduces the known accuracy for monomers, but deteriorates significantly 
with increasing chain length. On the other hand, the first two approaches provide rather good agreement. Particularly, 
the analytical result of Eq. (]38[) is surprisingly simple and yields predictions of similar good quality for all chain lengths 
studied. 

LGM wishes to thank Marcus Miiller, F. J. Bias, E. de Miguel, F. Bresme and C. Vega for helpful discussions; This 
work was supported by Ministerio de Educacion y Ciencia through a Ramon y Cajal grant and project FIS200766079- 
C02-00 and by Comunidad Autonoma de Madrid throught project MOSSNOHO-S0505/ESP/0299. PB acknowledges 
EU for partial funding this work as a TOK contract No. 509249. 



[1] J. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon, Oxford, 1982). 

[2] R. Seemann, M. Brinkmann, E. J. Kramer, F. F. Lange, and R. Lipowsky, Proc. Nat. Acad. Sci. 102, 1848 (2005). 

[3] A. Milchev, M. Muller, and K. Binder, Phys. Rev. E 72, 031603 (2005). 

[4] L. G. MacDowell and M. Muller, J. Chem. Phys. 124, 084907 (2006). 



13 



[5] J. S. Rowlinson, J. Stat. Phys. 20, 197 (1979). 

[6] J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 (1958). 

[7] P. Tarazona, Mol. Phys. 52, 81 (1984). 

[8] C. E. Woodward, J. Chem. Phys. 94, 3183 (1990). 

[9] M. Miiller and L. G. MacDowell, Macromolecules 33, 3902 (2000). 

[10] P. Bryk and S. Sokolowsky, J. Chem. Phys. 121, 11314 (2004). 

[11] P. Bryk, K. Bucior, S. Sokolowski, and G. Zukocinski, J. Physics-condensed Matter 16, 8861 (2004). 

[12] J. H. Irving and J. G. Kirkwood, J. Chem. Phys. 18, 17 (1950). 

[13] J. J. Magda, M. Tirrell, and H. T. Davis, J. Chem. Phys. 83, 1888 (1985). 

[14] A. Trokhymchuk and J. Alejandre, J. Chem. Phys. Ill, 8510 (1999). 

[15] P. Orea, Y. Duda, and J. Alejandre, J. Chem. Phys. 118, 5635 (2003). 

[16] E. M. del Rio and E. de Miguel, Phys. Rev. E 55, 2916 (1997). 

[17] F. Varnik, J. Baschnagel, and K. Binder, J. Chem. Phys. 113, 4444 (2000). 

[18] A. Milchev and K. Binder, J. Chem. Phys. 114, 8610 (2001). 

[19] D. Duque, J. C. Pamies, and L. F. Vega, J. Chem. Phys. 121, 11395 (2004). 

[20] K. Binder, Phys. Rev. A 25, 1699 (1982). 

[21] J. J. Potoff and A. Z. Panagiotopoulos, J. Chem. Phys. 112, 6411 (2000). 

[22] J. R. Errington, Phys. Rev. E 67, 012102 (2003). 

[23] L. G. MacDowell, J. Chem. Phys. 119, 453 (2003). 

[24] J. B. Hooper, J. D. McCoy, J. G. Curro, and F. van Swol, J. Chem. Phys. 113, 2021 (2000). 

[25] H. Furukawa and K. Binder, Phys. Rev. A 26, 556 (1982). 

[26] L. G. MacDowell, P. Virnau, M. Miiller, and K. Binder, J. Chem. Phys. 120, 5293 (2004). 

[27] M. Miiller and M. Schick, J. Chem. Phys. 105, 8282 (1996). 

[28] E. Chacon and P. Tarazona, J. Phys.: Condens. Matter 17, S3493 (2005). 

[29] H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 (1959). 

[30] J. R. Henderson, Mol. Phys. 50, 741 (1983). 

[31] P. Bryk, R. Roth, K. R. Mecke, and S. Dietrich, Phys. Rev. E 68, 031602 (2003). 

[32] F. Bresme and N. Quirke, Phys. Rev. Lett. 80, 3791 (1998). 

[33] F. Bresme and N. Quirke, J. Chem. Phys. 110, 3536 (1999). 

[34] G. J. Gloor, G. Jackson, F. J. Bias, and E. de Miguel, J. Chem. Phys. 123, 134703 (2005). 

[35] E. de Miguel and G. Jackson, Mol. Phys. 104, 3717 (2006). 

[36] H. L. Vortler and W. R. Smith, J. Chem. Phys. 112, 5168 (2000). 

[37] B. Widom, J. Chem. Phys. 39, 2808 (1963). 

[38] D. Chandler, J. D. McCoy, and S. J. Singer, J. Chem. Phys. 85, 5971 (1986). 

[39] E. Helfand, J. Chem. Phys. 56, 3592 (1972). 

[40] M. Miiller and L. G. MacDowell, J. Phys.: Condens. Matter 15, R609 (2003). 

[41] A. L. Frischknecht, J. D. Weinhold, A. G. Salinger, J. G. Curro, L. J. D. Frink, and J. D. McCoy, J. Chem. Phys. 117, 
10385 (2002). 

[42] J. B. Hooper, J. D. McCoy, and J. G. Curro, J. Chem. Phys. 112, 3090 (2000). 

[43] Y. X. Yu and J. Z. Wu, J. Chem. Phys. 117, 2368 (2002). 

[44] E. Kierlik and M. L. Rosinberg, J. Chem. Phys. 99, 3950 (1993). 

[45] M. S. Wertheim, J. Chem. Phys. 87, 7323 (1987). 

[46] R. Evans and U. M. B. Marconi, J. Chem. Phys. 87, 7138 (1987). 

[47] J. E. Finn and P. A. Monson, Phys. Rev. A 39, 6402 (1989). 

[48] J. E. Finn and P. A. Monson, Mol. Phys. 65, 1345 (1988). 

[49] J. Forsman and C. E. Woodward, Mol. Phys. 90, 637 (1997). 

[50] N. G. Almarza, E. de Miguel, and G. Jackson, private communication. 

[51] F. G. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001). 

[52] L. G. MacDowell, C. Vega, and E. Sanz, J. Chem. Phys. 115, 6220 (2001). 

[53] J. I. Siepmann and D. Frenkel, Mol. Phys. 75, 59 (1992). 

[54] L.-J. Chen, J. Chem. Phys. 103, 10214 (1995). 

[55] P. Orea, Duda, and Alejandre, J. Chem. Phys. 123, 114702 (2005). 

[56] J. K. Singh, D. A. Koike, and J. R. Errington, J. Chem. Phys. 119, 3405 (2003). 

[57] R. Dickman and C. K. Hall, J. Chem. Phys. 89, 3168 (1988). 

[58] A. Yethiraj and C. E. Woodward, J. Chem. Phys. 102, 5499 (1995). 

[59] M. Heni and H. Lowen, Phys. Rev. E 60, 7057 (1999). 

[60] F. Bresme, J. Phys. Chem. B 106, 7852 (2002). 

[61] H. Reiss, H. L. F. E. Helfand, and J. L. Lebowitz, J. Chem. Phys. 32, 119 (1960). 

[62] T. Boublik, J. Chem. Phys. 53, 471 (1970). 

[63] Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989). 

[64] R. Roth, R. Evans, A. Lang, and G. Kahl, J. Phys.: Condens. Matter 14, 12063 (2002). 

[65] Y. X. Yu and J. Z. Wu, J. Chem. Phys. 117, 10156 (2002). 

[66] P. M. Konig, R. Roth, and K. R. Mecke, Phys. Rev. Lett. 93, 160601 (2004). 

[67] H. Hansen-Goos and R. Roth, J. Phys.: Condens. Matter 18, 8413 (2006). 



14 



[68] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992). 

[69] P. Virnau and M. Miiller, J. Chem. Phys. 120, 10925 (2004). 

[70] M. Fitzgerald, R. R. Picard, and R. N. Silver, Europhysics Lett. 46, 282 (1999). 

[71] B. Smit, Mol. Phys. 85, 153 (1995). 



15 



TABLE I: Surface tension for the Spherically Truncated Lennard-Jones potential as determined by the wandering interface 
method with balance and bracket sampling. Results obtained from the virial method are also shown for comparison [14|] . Range 
column includes (ik for the counter tension method or the surface range sampled, for the bracket method. Size column indicates 
either the average lateral size for the wandering interface methods or the fixed lateral size for the virial method. Averages for 
balance sampling were collected over 400 Kcycles, except for those indicated with an asterisk, collected over twice as many. 
Averages for the bracket sampling were collected over 1000 Kcycles. 



k B T/e 


method 


constrain 


system size 




0.70 


balance 


0.276 


13.88 


0.781(8) 


0.70 


bracket 


[169, 181] 


13.19 


0.815(2) 


0.70 


bracket 


[179, 183] 


13.40 


0.797(10) 


0.70 


virial 


- 


13.41 


0.806(16) 


0.72 


balance 


0.175 


12.59 


0.75(2) 


0.72 


bracket 


[146, 158] 


12.11 


0.735(4) 


0.72 


virial 


- 


13.41 


0.743(3) 


0.80* 


balance 


0.200 


14.16 


0.594(5) 


0.80* 


balance 


0.150 


13.17 


0.598(6) 


0.80* 


balance 


0.125 


12.57 


0.601(6) 


0.80 


balance 


0.100 


11.83 


0.61(2) 


0.80 


bracket 


[179, 183] 


13.40 


0.591(7) 


0.80 


virial 


- 


13.41 


0.618(9) 


0.90* 


balance 


0.0756 


12.54 


0.413(3) 


0.90 


bracket 


[178, 184] 


13.38 


0.416(9) 


0.90 


virial 


- 


13.41 


0.438(4) 


0.92 


balance 


0.0903 


13.49 


0.377(7) 






n 77 18^1 


A-O.OO 


u.o i u yoj 


0.92 


virial 




13.41 


0.374(6) 


1.00 


balance 


0.0466 


13.14 


0.236(6) 


1.00 


bracket 


[172, 190] 


13.19 


0.246(5) 


1.00 


virial 




13.41 


0.244(5) 


1.127* 


bracket 


[133, 139] 


11.65 


0.0637(8) 


1.127 


bracket 


[170, 190] 


13.29 


0.0633(2) 


1.127* 


balance 


0.0057 


10.97 


0.070(3) 


1.127* 


balance 


0.0075 


11.76 


0.069(4) 


1.127 


virial 




13.41 


0.049(6) 
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TABLE II: Surface tension for the square well fluid with A = 1.5a. Results from this work are compared with the virial method 
(RefQJ) and with the Test Area Method (Ref|H). 



U Tic 

KB-L ft 


method 


constrain 


system size 


2 


n on 
u.yu 


bracket 


n iq 7p; 1 1 a op;1 
[1 (6. i O, 1 f 4.Z0J 


1 Q 1 o 


U.4ol(4UJ 


n on 

u.yu 


virial 




i n 
1U 


n /1oq/io\ 
U.42o(lU J 


n on 

u.yu 


1 AM 




i n ic< 

1U. ( O 


U.40o(lo ) 


n op; 

u.yo 


bracket 


n ic\ 7p; 1 71 op;1 

[1 f U. /D, 1 f l.ZOj 


1 Q HQ 


U.oDl(oo J 


0.95 


virial 




10 


0.349(9) 


n op; 
u.yo 


T'A 1\/T 
±AM 




1 n 52Q 

lU.oo 


U.ooO(zU j 


l.UU 


bracket 


ri 7fl 7 1 71 nl 

[1 /U. ( , 1 ( l.UJ 


1 Q H7 

lo.U f 


U.2DO(2U J 


1 nn 
l.UU 


virial 




i n 
1U 


U.zoo(lo ) 


i nn 
l.UU 


1AM 




1 1 m 
11. Ul 


U.2oo(2U J 


1.05 


bracket 


[170.15, 170.45] 


13.05 


0.166(21) 


1.05 


virial 




10 


0.202(11) 


1.05 


TAM 




11.18 


0.205(14) 


1.10 


bracket 


[168.75, 169.25] 


13.00 


0.132(22) 


1.10 


virial 




10 


0.142(12) 


1.10 


TAM 




11.39 


0.153(14) 



TABLE III: Wall-fluid interface tension for a bead-spring polymer model (FENE-LJ) at T=1.68. Interface tensions are given 
in units of the Lennard- Jones energy parameter e. 



C w 


"fwl 




3.30 


-0.181(6) 


-0.0139(3) 


3.20 


-0.164(12) 


-0.0119(3) 


3.10 


-0.136(15) 


-0.0106(3) 


3.00 


-0.120(10) 


-0.0095(6) 


2.80 


-0.098(7) 


-0.0072(1) 


2.00 


0.032(8) 


-0.0015(1) 


1.00 


0.161(10) 


0.0011(2) 


0.00 




-0.00009(10) 



TABLE IV: Chemical potentials (/i/fcsT) of tangent hard sphere chains for several chain lengths, m and monomer densities 
p. The results for m=12 and m=16 at p = 0.60(7~' j could not be measured from the NVT+test particle method and merely 
correspond to the chemical potential imposed during the simulations. Results for n/ksT are given as an excess over ideal 
chains. The entry for zero density is the difference between real and ideal chain contributions (see Refl7ll for a discussion on 
reference states appropriate to configurational bias grand canonical simulations). 



p/m 


2 


4 


8 


12 


16 


0.00 





0.6277 


2.1085 


3.6495 


5.2118 


0.10 


-2.34814(4) 


-2.12764(7) 


-0.8409(1) 


0.7589(2) 


2.4841(2) 


0.20 


-0.81907(5) 


-0.1513(1) 


1.9751(2) 


4.3995(6) 


6.9482(8) 


0.30 


0.67658(5) 


2.0178(3) 


5.4555(4) 


9.19(2) 


13.060(2) 


0.40 


2.4000(2) 


4.7275(3) 


10.123(2) 


15.81(7) 


21.59(7) 


0.50 


4.5358(6) 


8.296(2) 


16.5(1) 


25.19(3) 


34.(2) 


0.60 


7.308(1) 


13.126(5) 


25.2(4) 


43. 


54. 
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TABLE V: Wall-fluid interface tension for tangent hard sphere chains on athermal substrates for several chain lengths, m and 
monomer densities, p (cf. Tab II VI for chemical potentials). The entries of the table are interface tensions j/ksT. Results are 
given in units of the hard sphere diameter. 



p 


2 


4 


8 


12 


16 


0.10 


0.0447(5) 


0.0352(4) 


0.0295(4) 


0.0273(4) 


0.0248(2) 


0.20 


0.104(2) 


0.0905(10) 


0.0788(14) 


0.0755(10) 


0.0745(14) 


0.30 


0.197(6) 


0.159(6) 


0.148(4) 


0.149(5) 


0.153(3) 


0.40 


0.347(15) 


0.278(11) 


0.26(9) 


0.240(7) 


0.25(2) 


0.50 


0.44(3) 


0.43(2) 


0.40(3) 


0.375(11) 


0.40(2) 


0.60 


0.83(8) 


0.63(3) 


0.60(4) 


0.61(5) 


0.52(4) 
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FIG. 1: Liquid- Vapor surface tension of the spherically truncated (R c = 2.5) Lennard Jones model 7* = ji v a 2 /e as a function 
of temperature T* — ksT/e. The squares and circles are obtained from the contertension and the bracketed sampling methods, 
respectively. The stars are results from |l4( |. The inset is a plot of inverse squared surface area against l//3£, with slope equal 
to f3j. 
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FIG. 2: Liquid- Vapor surface tension of a square well fluid with A = 1.5c. Results are given in reduced units, with 7* = 7;„cr 2 /e 
plotted as a function of temperature T* = ksT/e. Circles, results from the bracketed sampling method of this work. Squares, 
anisotropy of the virial tensor (l5| . Diamonds, Test Area Method [Sj. Triangles, Binder's method [5r3 |. 
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FIG. 3: Wall-fluid interface tensions for the LJ-FENE model. The figure shows the spreading coefficient ^ wv —"/ w i as a function 
of wall strength, e m . Squares are results obtained in this work from explicit evaluation of y wv and y w i. Circles and lines are 
results for the spreading coefficients from Ref[^ The dashed horizontal lines correspond ±7;„ and the intersections correspond 
to wetting and drying transitions. The inset shows the tensions as a function of wall strength. Circles and squares refer to the 
wall-liquid and wall- vapor tensions, respectively. 
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FIG. 4: Surface area mean squared displacement as a function of Monte Carlo cycles for chains of length m=16. From top 
to bottom, results refer to monomer densities p = O.lOcr^ 3 (full line), p = 0.30cr -3 (dashed line, 2 times magnified) and 
p — 0.50cr -3 (dash-dotted line, 10 times magnified). 
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FIG. 6: Wall-fluid interface tensions for hard sphere chains of length m = 12 adsorbed on an athermal substrate with the 
dividing surface located at z = \a. The circles are simulation results with the volume of the slit pore defined as V = (L z — a). 
The filled squares are results obtained by substracting |pcr to the simulated data for 7 z =o- The lines are theoretical results. 
Full line, TPT1-DFT; dashed line, SPT; dotted lines, DFT-based analytical approximation. 



